# !/usr/bin/env python
# -*- coding: utf-8 -*-
"""
@Time        : 2020/11/27 12:12
@Author      : Albert Darren
@Contact     : 2563491540@qq.com
@File        : integrate.py
@Version     : Version 1.0.0
@Description : TODO
@Created By  : PyCharm
"""
from scipy.integrate import simps, trapz, romberg
import numpy as np
from sympy.abc import x
from sympy import exp


# 被积函数1
def f1(x_coordinates):
    # 此处会有运行时警告，因为自然对数0处无意义
    y = np.log(x_coordinates) * np.power(x_coordinates, 0.5)
    # 积分结点0函数值为nan无意义，需要用理论特殊处理
    if x_coordinates[0] == 0:
        y[0] = 0
    return y


# 被积函数2
def f2(x_coordinates):
    # 此处会有运行时警告，因为自然对数0处无意义
    y = np.log(x_coordinates) * np.power(x_coordinates, 0.5)
    return y


# 积分精确值
exact = (-4 / 9)
print("精确值为:{}".format(exact))
for n in [10, 20]:
    # 积分结点横坐标和步长
    x, step = np.linspace(0, 1, n + 1, retstep=True)
    # 复合梯形求积公式
    com_trapz_integral = trapz(f1(x), x)
    error = abs(exact - com_trapz_integral)
    print("""步长h={},复合梯形求积公式积分为:{},误差为:{}""".format(step, com_trapz_integral, error))
    # 复合辛普森求积公式
    com_simpson_integral = simps(f1(x), x, even='avg')
    error = abs(exact - com_simpson_integral)
    print("""步长h={},复合辛普森求积公式积分为:{},误差为:{}""".format(step, com_simpson_integral, error))
# 龙贝格求积公式,由于被积函数在0处无意义，无法直接使用0作为积分下限，只能用任意极小值逼近0
romberg_integral = romberg(f2, 1e-20, 1, divmax=20)
error = abs(exact - romberg_integral)
print("""外推次数20次，龙贝格求积公式积分为:{},误差为:{}""".format(romberg_integral, error))
